
install.packages("car")

library(car)


#setwd("~/Dropbox/Job Market Paper/IO Revisions/Replication_Files/Replication_Appendix/A_Fig_1_2")


data_q<-read.csv(file="data_fig_A1_A2.csv")


###################### Part 1 of Appendix ######################

###################### Figures # A2 #### ######################

pdf(file="AFigure_2.pdf")


eleventhun<-subset(data_q,data_q$Year==1100)

eleventhun50<-subset(data_q,data_q$Year==1150)

twelvehun<-subset(data_q,data_q$Year==1200)

twelvehun50<-subset(data_q,data_q$Year==1250)

thirteenhun<-subset(data_q,data_q$Year==1300)

thirteenhun50<-subset(data_q,data_q$Year==1350)

fourteenhun<-subset(data_q,data_q$Year==1400)

fourteenhun50<-subset(data_q,data_q$Year==1450)

fifteenhun<-subset(data_q,data_q$Year==1500)

fifteenhun50<-subset(data_q,data_q$Year==1550)

sixteenhun<-subset(data_q,data_q$Year==1600)

sixteenhun50<-subset(data_q,data_q$Year==1650)

seventeenhun<-subset(data_q,data_q$Year==1700)

seventeenhun50<-subset(data_q,data_q$Year==1750)

seventeenhun90<-subset(data_q,data_q$Year==1790)

par(mfrow=c(2,7),mar=c(1,1,1,1),oma=c(1,1,1,1),omi=c(1,1,1,1))


qqnorm(scale(twelvehun$Area),main="1200",xlab="",ylab="Sample Quantile of Area",)
abline(0,1)



qqnorm(scale(thirteenhun$Area),main="1300",xlab="",ylab="")
abline(0,1)


qqnorm(scale(fourteenhun$Area),main="1400",xlab="",ylab="")
abline(0,1)


qqnorm(scale(fifteenhun$Area),main="1500",xlab="",ylab="")
abline(0,1)


qqnorm(scale(sixteenhun$Area),main="1600",ylab="",xlab="")
abline(0,1)


qqnorm(scale(seventeenhun$Area),main="1700",ylab="",xlab="")
abline(0,1)


qqnorm(scale(seventeenhun90$Area),main="1790",ylab="",xlab="")
abline(0,1)

##


qqnorm(scale(log(twelvehun$Area)),main="",xlab="",ylab="Sample Quantile of log(Area)",)
abline(0,1)



qqnorm(scale(log(thirteenhun$Area)),main="",xlab="",ylab="")
abline(0,1)


qqnorm(scale(log(fourteenhun$Area)),main="",xlab="",ylab="")
abline(0,1)


qqnorm(scale(log(fifteenhun$Area)),main="",xlab="",ylab="")
abline(0,1)


qqnorm(scale(log(sixteenhun$Area)),main="",ylab="",xlab="")
abline(0,1)


qqnorm(scale(log(seventeenhun$Area)),main="",ylab="",xlab="")
abline(0,1)


qqnorm(scale(log(seventeenhun90$Area)),main="",ylab="",xlab="")
abline(0,1)



mtext("Normal Quantiles",side=1, outer=TRUE, line=4) 

mtext("Sample Quantiles",side=2, outer=TRUE, line=4) 

mtext("QQ Plots",side=3, outer=TRUE, line=4,cex=2) 


dev.off()

###################### Figures # A1 #### ######################

pdf(file="AFigure_1.pdf")

par(mfrow=c(2,7),mar=c(1,1,1,1),oma=c(1,1,1,1),omi=c(1,1,1,1))


density_a<-density(twelvehun$Area)


plot(density_a,xlab="",ylab="Density",axes=F,main="1200")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)


mtext(expression(KM^2),2,line=3)



density_a<-density(thirteenhun$Area)


plot(density_a,xlab="",ylab="Density",axes=F,main="1300")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)




density_a<-density(fourteenhun$Area)


plot(density_a,xlab="",ylab="Density",axes=F,main="1400")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)





density_a<-density(fifteenhun$Area)


plot(density_a,xlab="",ylab="Density",axes=F,main="1500")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)







density_a<-density(sixteenhun$Area)


plot(density_a,xlab="",ylab="Density",axes=F,main="1600")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)






density_a<-density(seventeenhun$Area)


plot(density_a,xlab="",ylab="Density",axes=F,main="1700")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)





density_a<-density(seventeenhun90$Area)


plot(density_a,xlab="",ylab="Density",axes=F,main="1790")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)






density_a<-density(log(twelvehun$Area))


plot(density_a,xlab="",ylab="Density",axes=F,main="")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)


mtext(expression(log(KM^2)),2,line=3)



density_a<-density(log(thirteenhun$Area))


plot(density_a,xlab="",ylab="Density",axes=F,main="")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)




density_a<-density(log(fourteenhun$Area))


plot(density_a,xlab="",ylab="Density",axes=F,main="")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)





density_a<-density(log(fifteenhun$Area))


plot(density_a,xlab="",ylab="Density",axes=F,main="")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)







density_a<-density(log(sixteenhun$Area))


plot(density_a,xlab="",ylab="Density",axes=F,main="")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)






density_a<-density(log(seventeenhun$Area))


plot(density_a,xlab="",ylab="Density",axes=F,main="")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)





density_a<-density(log(seventeenhun90$Area))


plot(density_a,xlab="",ylab="Density",axes=F,main="")

axis(1); axis(2)

polygon(c(0,density_a$x),c(0,density_a$y),col="black",density=3,angle=45,border="black",fillOddEven=T)




mtext("State Size Kernel Density",3,outer=T,line=3,cex=2)


dev.off()


############## Box-Cox Transform ##################


print("This gives the box-cox transformation from the appendix")

lam<-powerTransform(data_q$Area)
lam<-coef(lam)


### This gives the result rho = -.099 ###

print(lam)

